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§ ■ ABSTRACT 
(N 

Global axisymmetric stability of viscous, resistive, magnetized Couette flow is re- 
' examined, with the emphasis on flows that would be hydrodynamically stable according 

to Rayleigh's criterion: opposing gradients of angular velocity and specific angular 
momentum. In this regime, magnetorotational instabilities [MRI] may occur. Previous 
work has focused on the Rayleigh-unstable regime. To prepare for an experimental study 
of MRI, which are of intense astrophysical interest, we solve for global linear modes in 
a wide gap with realistic dissipation coefficients. Exchange of stability appears to occur 
through marginal modes. Velocity eigenfunctions of marginal modes are nearly singular 
£^ ' at conducting boundaries, but magnetic eigenfunctions are smooth and obey a fourth- 

order differential equation in the inviscid limit. The viscous marginal system is of tenth 
order; an eighth-order approximation previously used for Rayleigh-unstable modes does 
not permit MRI. Peak growth rates are insensitive to boundary conditions. They are 
q . predicted with surprising accuracy by WKB methods even for the largest-scale mode. 

We conclude that MRI is achievable under plausible experimental conditions using easy- 
^ | to-handle liquid metals such as gallium. 

> 

X. 

1. Introduction 

To an even greater extent than large-scale terrestrial ones, astrophysical flows are nearly invis- 
cid. Yet observations show that they dissipate efficiently. For example, accretion disks (flattened 
systems of gas in orbit about a star or black hole) must lose orbital energy in order that the gas flow 
onto the central object. The influence of turbulence has long been suspected, but purely hydrody- 
namic turbulence is probably ineffective because of the strongly stable radial angular momentum 
gradient in the disk. When warm enough to be partially ionised, as they often are, accretion disks 
become magnetohydrodynamic (MHD) fluids. It is now believed that turbulence and orbital de- 
cay are driven by magnetorotational instabilities (MRI). Although discovered by Velikhov (1959) 
and Chandrasekhar (1960), MRI did not come to the attention of the astrophysical community 
until rediscovered by Balbus & Hawley (1991). MRI requires that the angular velocity (O) must 
decrease with distance from the axis, but it is distinguished from Rayleigh's centrifugal instability 
by persisting when the gradient of specific angular momentum is positive, i.e. <9(r 4 f2 2 ) / dr > 0. A 
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weak background magnetic field is required, and MRI can be axisymmetric when the field parallel 
to the rotation axis. The rotation law of accretion disks satisfies both conditions above (usually 
oc r 1 / 2 ), and while the magnetic geometry is uncertain, it has little influence on the maximum 
MRI growth rate in ideal MHD unless the field is purely toroidal (Balbus & Hawley 1992; Terquem 
& Papaloizou 1996). 

There exists a body of experimental work on magnetized Couette flow (Donnelly & Ozima 1960, 
1962; Donnelly & Caldwell 1964; Brahme 1970), but the MRI has never been demonstrated in the 
laboratory. The main obstacle is that liquid metals are strongly resistive on laboratory scales, with 
magnetic diffusivity r\ > 10 3 cm 2 s _1 . The viscosity is much smaller, typically v ~ 10 -3 cm 2 s _1 . 
Their ratio is the magnetic Prandtl number P m = u/r] ~ 10~ 6 . 

Chandrasekhar (1961) also analyzed dissipative magnetized Couette flow. After laying out 
general equations for viscous and resistive linearized axisymmetric perturbations, he discarded a 
term involving shear from the azimuthal induction equation on the grounds that P m <C 1. In 
fact the neglected term involves neither v nor rj directly. Thus the commonly used name "small- 
P m approximation" is somewhat unfortunate. It presumes that viscous and inertial forces are 
comparable throughout the flow. Although this is often the case for marginal Rayleigh instabilities 
resisted by viscous (and perhaps magnetic) forces, it is not the case for the MRI modes of interest 
here, where viscosity is important in boundary layers only. We shall show that MRI modes do not 
exist in Chandrasekhar's approximation (§3). 1 We shall introduce a different P m — > limit that 
retains all terms in the induction equation but neglects viscous boundary layers. 

Chandrasekhar's results have been refined by Chang & Sartory (1967); Hassard, Chang, & Lud- 
ford (1972); Vislovich, Novikov, & Sinitsyn (1986); Takhar, Ali, & Soundalgekar (1989); Soundal- 
gekar, Ali, h Takhar (1994); Chen & Chang (1998), but always under Chandrasekhar's "small-P m " 
approximation. 

In a previous paper (Ji, Goodman, & Kageyama 2001, henceforth Paper I), we have used local 
WKB methods to survey the MRI regime for realistic materials and laboratory parameters. The 
most unstable modes (and perhaps the only ones accessible to experiment in the near term) have 
wavelengths twice as large as the apparatus, so that WKB methods are not to be trusted a priori. 
The present paper therefore discusses the global linear analysis. We integrate the full set of viscous 
and resistive equations via an initial-value scheme to obtain numerical growth rates for cases that 
would be stable by Rayleigh's criterion. Boundary conditions are problematic for inviscid marginal 
modes. Nevertheless, the locally obtained growth rates are found to be good approximations, even 
though the radial eigenfunctions are far from being sinusoidal. We predict instability under feasible 
conditions: gap widths and heights of order ten centimeters, fields strengths of several kilogauss, 
and rotation rates of several hundred radians per second. 

A supplementary analysis shows that curves of marginal MRI stability are well approximated 



x Of course he did find MRI modes, but in a separate analysis assuming ideal MHD (Chandrasekhar 1960). 
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by completely neglecting viscous terms, especially for insulating cylinders. Although P m — > 0, this 
is not Chandrasekhar's approximation. It eliminates the velocity perturbations from the analysis 
and leads to a system of only four radial derivatives. The boundary conditions on velocity are 
not satisfied by the reduced system even with stress- free (slipping) boundaries. This is resolved by 
restoring the viscosity but there is little change in the magnetic eigenfunctions or the growth rate 
when P m is realistically small. Hence when MRI are of primary interest, dimensionless numbers 
based on viscosity, viz. Hartmann and Taylor numbers, are less useful for characterizing stability 
than numbers that remain finite as P m — ► 0, such as Lundquist number and magnetic Reynolds 
number. 



2. Basic equations and boundary conditions 

We use cylindrical coordinates r9z aligned with the rotation of the fluid, and gaussian elec- 
tromagnetic units. In equilibrium, the magnetic field is constant and parallel to the axis; it is 
described by the associated Alfven speed Va = B /^4ttp, where p is the (constant) density of the 
liquid metal. The angular velocity is Cl(r). Perturbations are axisymmetric and sinusoidal in z 
with wavelength 2n/k: 

5B r /^4irp = p r (r, t) cos kz Sv r = ip r (r, t) sin kz 

5Bq I Airp = Pe(r,t) cos kz 5vg = ipg(r,t)smkz (1) 

5B z /y/4-Kp = (3 Z ( r ) t) sin kz 5v z = <p z (r, t) cos kz 

We often write h = ir/k for half the wavelength, having in mind an experiment of finite height h 
and rigid vertical boundaries. 

Perturbations associated with a single mode are exponential in time, but it is convenient 
to allow for general time dependence and discover the fastest-growing mode by integrating the 
linearized equations forward in time. These equations are (see the Appendix) 



06 = v(D-k 2 )p e + kV A iP8 + rn% (2) 

<p e = u(D - k 2 )ip e - WaPo - r -1 (r 2 fi)Vr (3) 

P r = f](D-k 2 )p r + kV A ip r (4) 

fa = v{D - k 2 )ip r - kV A p r + n (5) 

(k 2 -D)U = 2Qk 2 ip e (6) 



The prime in eqs. (2)-(3) means d/dr, and 

D = — -—-1 

dr 2 r dr r 



Note that Chandrasekhar (1961) denotes the latter operator by DD*. 
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Equations (2)-(3) are the azimuthal components of the induction and euler equations, while 
(4)-(5) are the corresponding radial components. By eliminating the auxiliary function II between 
eqs. (5)-(6), one obtains an equation equivalent to Chandrasekhar's eq. (168), although some dif- 
ferences in sign occur because we have taken perturbations proportional to sin kz where he took 
coskz and vice versa. Our eqs. (2), (3), and (4) correspond to his eqs. (163), (160), and (162), 
respectively. 

The flow is confined between concentric cylinders with radii r± < r<i- If these are perfectly 
conducting, the magnetic boundary conditions are 



(3 r = 0, {r/3 e y = 0. 

If perfectly insulating, then (I n and K n are modified Bessel functions) 



(7) 



(3 r x < 
Pe = 0. 



[krIo(kr)]/Ii(kr) at r = r\ 
— [krKo(kr)]/ K±(kr) at r = T2 



The conditions on velocity are 



(fir = 0, 



vipo = = v(r(p r )'. 



(8) 



(9) 



We have put the viscosity in the latter two conditions so that tpg and ip z = (rip r )'/(kr) will be 
unconstrained when the flow is inviscid. 

The insulating conditions (8) are not accurate for an experiment of finite height h = vr/fc, 
since they assume a vertically periodic solution for the vacuum field outside the cylinders. The 
conducting conditions (7) do not have this drawback. In both cases, there will be viscous boundary 
layers at the top and bottom (unless the end caps rotate differentially), but we expect that the 
error committed by neglecting those layers is small for P m ~ 10 -6 and growth times shorter than 
the Ekman circulation time. In any event, the end caps should be insulating so that no magnetic 
stress acts upon them. 



3. Why the small-P m approximation suppresses the MRI 

The underlined term in eq. (2) is the critical one that Chandrasekhar (1961) and subsequent 
authors neglected on the grounds that P m <C 1. To see that this term is necessary to the MRI, it 
is useful to reduce eqs. (2)-(6) to a single equation. 

For brevity's sake, assume a mode with definite growth rate s, and write = D — k 2 , 
loa = kVA- Eq. (4) yields 

= ^a 1 (s - f]D k )f3 r . (10) 
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Substituting for <p r in eq. (5), applying D kl and eliminating D^H via eq. (6), one has 

<pe = -(2Qk 2 LO A y 1 [(a - vD k ){s - r)D k ) + lo 2 a ] D k f3 r . 
Solving for (3g from eq. (3) and eliminating ipg and ip r in favor of (3 r , one has 

Po = |(s - "AO^fcS [( s ~ ^fc)( s - V D k) +^a] D k 

Using these to eliminate <pg and fig from eq. (2) yields the desired tenth-order equation: 
| [(s - uD k ){s - r)D k ) + u 2 A ] ^ [(a - i/D fc )(a - 7/D fc ) + u>\] (-k' 2 D k ) 



+ (s-r ] D k )^^( S -7 ] D k )\(3 r = -io 2 A rn'(3 r . 

r ' - 



(11) 



(12) 



(13) 



We are interested in Rayleigh-stable cases. Without loss of generality, we may assume that the 
angular velocity (Q) and vorticity [r -1 (r 2 Sl)'] are positive throughout the flow, but the shear (rW) 
may be negative. Now D k is clearly negative definite and self adjoint with either of the boundary 
conditions (7) or (8). In the narrow-gap limit (r2 — r±)/(r2 + r±) — > 0, the angular velocity and 
vorticity are positive constants. It follows that the operator in braces in the lefthand side of eq. (13) 
is positive definite for nonnegative real s. Therefore, at least in the narrow-gap limit, there can be 
no modes with positive real growth rate when the underlined shear term is neglected and the Rayleigh 
stability criterion is satisfied. We interpret this to mean that the MRI is not present. 

Previous studies of the time-dependent problem have neglected the time derivatives in the 
induction equation as well as the underlined term (Chang & Sartory 1967; Chen & Chang 1998); 
in this case, the operator in question becomes quadratic in s so that the coefficient of Imag(s) is 
positive-definite for Real(s) > 0, which rules out complex growing modes, i.e. overstabilities. 

Unfortunately, we cannot draw such strong conclusions in the wide-gap case where Q, and 
perhaps also r _1 (r 2 il)' vary significantly with radius. The operator on the left side of eq. (13) is no 
longer self adjoint in general, because D k and Q do not commute. For marginal modes, however, 



uDi + 



2 1 



2Qk 2 



uDi + 



2 1 



V D k 



{r 2 ny 



<Pr 



(14) 



in which we have used eq. (10). In this case the problem is only eighth order in radial derivatives if 
the righthand side is neglected, as noted by Chandrasekhar (1961). More importantly, the operator 
in these curly braces is self adjoint in the interesting special case that 



Q(r) 



a + 



,2 ' 



(15) 
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since r~ l (r 2 Q)' = 2a is then constant, and the rest of the operator is symmetrical. Eq. (15) is 
the angular-velocity profile of a Couette flow in steady state, because it implies a radially constant 
viscous angular-momentum flux. We conclude that there are no Rayleigh-stable marginal modes 
when the magnetic shear term is neglected, even for wide gaps. 



4. Marginal inviscid modes 

If v = 0, then as the growth rate s — > 0, eq. (13) reduces to 



'2 D (r 2 ny u>\ 



D k p r = -u\rQ!_§r. (16) 



Even with the righthand side included, this is only a fourth-order system. Paradoxically, there are 
six boundary conditions to be satisfied: (pg and ip z are not constrained when v = 0, but ip r = at 
both boundaries, in addition to a total of four magnetic conditions. For the insulating case (8), the 
paradox is resolved because eqs. (10) and (12) show that both ip r = and (3g = are equivalent 
to Dj,(3 r = 0, so that there are only four independent boundary conditions after all. But in the 
conducting case (7), we have via eqs. (10)-(12) that (3 r = Dk(3 r = (rDkj3 r )' = at both boundaries, 
and not all six conditions can be satisfied. 

The crux of the difficulty is the azimuthal euler equation (3), which reduces to an algebraic 
relation. At zero viscosity and growth rate, azimuthal force balance requires 

B d5Bg 

™V r = —— , (17) 

so that the Lorentz force (Sj r x Bq) balances the Coriolis force, which vanishes at the boundary. If 
the boundary is insulating, then SBg (oc Sj r ) also vanishes. But at conducting boundaries, 5Bg ^ 
(Sj r 0), so that viscosity must intervene to maintain azimuthal force balance. For small P m , the 
marginal eigenfunction displays a dramatic boundary layer (Fig. 3). 

Viscous boundary layers are common in hydrodynamics. Normally they occur because the 
tangential component of velocity must match that of the boundary itself, even when the viscosity is 
small (a "no-slip" boundary condition). For conducting cylinders, however, a boundary layer would 
occur even if the viscous stress vanished at the boundary, because of the impossibility of satisfying 
eq. (17). In the present case, the viscous layer is driven by tangential magnetic field (or normal 
component of current) rather than tangential velocity. 

To summarize, the inviscid limit is singular for marginally stable modes. The eigenfunctions 
become ill-behaved because there are more boundary conditions than radial derivatives to satisfy 
them, at least for conducting boundaries. The numerical evidence presented below indicates, how- 
ever, that the locus of marginal stability in the parameter space of equilibria is actually continuous 
as P m — > 0. Hence the relatively simple fourth-order differential equation (16) predicts marginal 
stability reasonably well when P m is sufficiently small. 



5. Numerical results 



We have approximated eqs (2)- (6) by finite-difference equations on a radial grid uniformly 
spaced in x = lnr. The background angular velocity has the form (15), since this is easiest to 
realize experimentally. The grid spacing Ax must be chosen fine enough so that Ar 2 /v is smaller 
than the interesting physical timescales in the problem, viz. u^ 1 and (r/d 2 )^ 1 (where d = r<i — t\ is 
the gap width), otherwise the viscous boundary layer will not be resolved. The minimum number 
of grid cells N ~ P m 1 ^ 2 ~ 10 3 . Our scheme has second-order spatial accuracy, even at the the 
boundaries. 

To ensure numerical stability, we use fully implicit time differencing. Spatial differences are 
written in terms of the unknown dependent variables at the new time step, so that a large linear 
system must be solved. Actually, our finite-difference matrix is band diagonal with 10 nonzero 
codiagonals in each of 5N rows, and it is independent of time step. We perform LU decomposition 
at the start of the evolution so that only the back substitution must be performed anew at each 
step. 

When a growing mode is present, it eventually dominates. Then the perturbation in radial 
magnetic field at successive time steps t n and t n+ \ = t n + At are related by, for example, 

(1 - sAt) p r ( Xj , t j+1 ) = 0r(Xj,tj), (18) 

if s > is the appropriate eigenvalue of the matrix defined by the spatial difference scheme and 
therefore an estimate of the physical growth rate, and (3 r is the corresponding eigenvector. Given 
$ r (xj,tj) and /3 r (xj,tj + i), we can compute s from eq. (18) without any truncation error in At as 
long as sAt < 1. The eigenfunctions are similarly independent of At. This is advantageous close 
to marginal stability where s is small. Our procedure is equivalent to finding the most-positive 
eigenvalue of the time-evolution matrix by inverse iteration. By extending eq. (18) to a three-term 
recurrence relation, we have allowed for complex eigenvalues. But in practice, all of our growing 
modes appear to have purely real values of s. Of course, our method is not immune to spatial 
truncation errors; these are 0(Ax 2 ) because we use second-order spatial differencing. 

Our initial-value code can treat slightly stable cases as well as unstable ones. By interpolation, 
we find parameters for marginal stability. Results are shown in Figs. (1) & (2) for material properties 
approximating liquid gallium, viz. p = 6g cm" 3 , rj = 2000cm 2 s~ 1 and P m = 1.6 x 1CT 6 . 

Marginal stability defines one constraint among the eight parameters defining the Couette flow: 
77, v, n, r2, k, S7i, Q2, and Va- Six independent dimensionless combinations of these can be formed. 
The magnetic Prandtl number P m = u/r] is one of these. The aspect ratio 

As!1 ir 1 - < 19 > 

and the elongation of the toroidal cross section, 
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define the geometry, where d = — r\ is the gap width and h = ir/k is half the vertical wavelength. 
It turns out that the magnetic eigenfunctions (3 r and (5q of the most unstable mode have a roughly 
parabolic dependence on r, and since one or the other vanishes at both boundaries, the effective 
horizontal wavenumber is ~ ir/d. The total wavenumber is then 



The Lundquist number 



K^le + ^f-lyfc?. (21) 

Wa = VAh J2_ 
~ V K 2 27T?? e 2 + 1 1 ) 

scales the Alfven frequency against the magnetic diffusion rate r]K 2 . In astrophysics, S is often 
called "magnetic Reynolds number." The plasma community generally reserves the latter term for 
a quantity involving fluid velocity, so we define the local magnetic Reynolds number by 

R m (r) = (23) 

The viscous Reynolds number is of course R m /P m . The dimensionless vorticity 

CM-*** 

parametrizes the angular momentum gradient, and the Rayleigh stability criterion is simply £(r) > 
0. In the astrophysical literature, the radial variation of angular velocity is often described by an 
index 

dlntt 

q = — so that Q = 2 — q. 

amr 

Of course, Q = 2 and q = in a uniformly rotating flow. 

When the aspect ratio is modest, R m and ( may vary considerably across the gap, and it is 
useful to define mean values of these dimensionless parameters. Following Paper I, we introduce 
Cl = and 

Cl - 2(rfo 2 -rfoi) 

Km -^> C - (r 2 2 -rf)0 (25j 

The locus of marginal stability is actually a hypersurface in the space (P m , A, e, £, S, Rm)- The 
curves in Figs. (l)-(2) are cuts through this locus at constant values of the first four parameters: 
P m = and P m = 1.6 x 10~ 6 ; A = e = 1; and two positive values of £ as indicated. The curves are 
drawn in physical units for the density and diffusivity of gallium. 

The inviscid results shown in in these figures were calculated by an independent numerical 
method based on eq. (16), which we have unpacked into a pair of second-order equations [using 
eq. (12) with s,v -> 0]: 

D <* = 5$y* = K ''wJ°- <26> 
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n R ^ r r rfl ' R K 2 

Dk(3e - 2WF A " ~ 



(l + e 2 )S 4 



■fa + (2 - QRmPr 



(27) 



2(R> 

together with the magnetic boundary conditions (7) or (8). Because of the large magnetic diffu- 
sivity, the magnetic variables are very well-behaved, so that this fourth-order system can be solved 
efficiently by a shooting method. 

Comparing Figs. 3 and 4, one sees that viscosity is more important for conducting boundary 
conditions than for insulating ones. In the conducting case, (p r and f}$ are nearly proportional to 
one another throughout most of the flow [Fig. 3]. This follows from eq. (3) in the limit s, v — > 0, as 
we discussed in §4. Since the radial velocity perturbation (oc ip r ) must vanish at the wall but the 
azimuthal magnetic perturbation (oc (3g) does not, there is a thin boundary layer in which viscous 
stress balances the azimuthal magnetic force. The righthand panel shows that the boundary layer 
is well resolved by these calculations, which use 4000 grid points uniformly spaced in lnr across 
the gap. At an insulating boundary, on the other hand, fig vanishes with ip r , and this leads to a 
much less dramatic viscous layer (Fig. 4). 

The narrow-gap limit A — > oo is experimentally impractical but theoretically important. Fig. 5 
shows eigenfunctions and curves of marginal stability and in this limit. Because of the boundary 
conditions, the eigenfunctions cannot be sinusoidal in r or x = (r — r±)/d, even though the equations 
of motion have elementary solutions of this form. The fourth-order inviscid system (26)- (27) has 
four roots for the radial wavenumber (or two if sign is ignored) at given parameters (e, S, R m , Q of 
the equilibrium, which are constant across the gap. (Note K is not an independent parameter, since 
Kd = \/e 2 + 1.) The solutions satisfying the boundary conditions are linear combinations of four 
complex exponentials, each containing one of these wavenumbers, and the parameters (e, S, R m , () 
must satisfy one constraint in order that a solution exist. If, however, one simply sets k r = ir/d in 
hopes of obtaining an approximate constraint, then — ► — K 2 and eqs. (26)- (27) or (16) would 
yield 

R ™ = 52 2(2-C-C5- 2 )' (28) 
This corresponds to the dispersion relation for marginal modes obtained in Paper I from a local 
WKB analysis with a ratio e of vertical to horizontal wavelength. Evidently there exists a minimum 
Lundquist number for instability, 

; s (29) 



and in the opposite limit of large S, 



i?m - Q - / 1 + 62 rso) 

Fig. (5) shows that the predictions (28)-(30) are qualitatively correct. 

Numerical growth rates are given in Table 1 for two representative angular-velocity profiles, 
both with T2 = Sr\ = 15 cm, h = 10 cm, and the material properties of gallium. The first case has 
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C = 0.063 and hence is Rayleigh-stable. The second has C = —0.019, so that the Rayleigh instability 
occurs at zero field strength. Growth rates are shown for conducting and insulating radial boundary 
conditions. Larger fields are required to initiate and to quench MRI with insulating boundaries 
than with conducting ones, presumably because perturbed lines of force expand past insulating 
walls into a volume slightly larger than the Couette flow itself. The final column of this table shows 
growth rates computed from the algebraic local dispersion relation given by in Paper I. Once again 
the WKB analysis predicts the global growth rates remarkably well, even though the wavelengths 
involved are actually larger than the gap width, and the angular velocity and shear rate vary by a 
factor « 9 across the gap. 

Figures (6) and (7) show two-dimensional cross sections of selected modes from Table 1. The 
flux and stream functions are related to the poloidal perturbations by 

i^L = V0xVx, 5v P = V0xVip. (31) 

In all cases, the poloidal velocity field consists of a single roll. The effect of the choice of boundary 
conditions is seen most clearly in the toroidal perturbations. In the first part of Fig. (6) and the 
second part of Fig. (7), it looks as though Svg does not vanish at the inner boundary; in fact it 
does vanish, but the viscous boundary layer is too small to be resolved by these plots. The first 
part of Fig. (7) does not show this behavior because the magnetic forces are absent; this mode is 
a classical hydromagnetic centrifugal instability. A boundary layer of the ordinary nonmagnetic 
variety occurs in Sv z . The corresponding magnetized case shown in the second part of the figure 
has a magnetically-driven boundary layer similar to that of the Rayleigh-stable flows. 



6. Summary and discussion 

We have presented a global linear stability analysis for magnetized Couette flow, including 
the dissipative effects of viscosity and resistivity, in regimes where magnetorotational instability 
(MRI) is possible. In view of the actual properties of liquid metals, and for plausible experimental 
lengthscales, resistivity is the main obstacle that must be overcome to demonstrate MRI. Previous 
theoretical studies of magnetized Couette flow have focused on the problem of suppressing Rayleigh 
instability with a magnetic field, and they have simplified the induction equations so as to reduce 
the number of radial and time derivatives in the problem. Such approximations may be adequate 
for the regime of interest to those studies, but they will not do if one is interested in MRI. Therefore 
we have worked with the full induction equations. 

Particular attention has been given to marginal stability. All of our numerical evidence indi- 
cates that exchange of stability occurs at vanishing complex growth rates; we have not encountered 
any overstable modes. In the inviscid limit of marginal stability, the number of radial derivatives 
reduces from ten to four, but the velocity perturbations become singular at conducting boundaries. 
The dominant singularity arises from an unbalanced tangential magnetic force, not from the usual 
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pressure and inertial effects that cause boundary layers in unmagnetized fluids. Despite these com- 
plications, the inviscid approximation predicts the locus of marginal stability reasonably well for 
liquid metals. 

Solving the linearized initial-value problem by finite-differences, we have calculated growth 
rates and stability boundaries for a liquid metal approximating gallium in an experimentally plau- 
sible geometry. For easier comparison with other theoretical work, we have also made calculations 
in the narrow-gap limit and expressed our results in the dimensionless coordinates of Lundquist 
number and magnetic Reynolds number. 

Remarkably, the growth rates are reasonably well predicted by a simple WKB approximation 
even though the WKB modes do not satisfy the boundary conditions and have a radial wavelength 
twice the gap width, and even when the rotation rate varies by an order of magnitude across the 
gap. We conclude from this that the algebraic WKB dispersion relation can be used for preliminary 
experimental design, at least for aspect ratios no more extreme than considered here (v2 ■ r\ = 3 : 1). 

There are good reasons to attempt an MRI experiment. First, one can hardly exaggerate the 
importance of this instability: few or no plausible alternative explanations exist for the dissipation 
of orbital energy in accretion disks, which are fundamental to so many of the most energetic sources 
known in the universe. Yet all present knowledge of the instability is purely theoretical, based as it is 
on linear analysis and computer simulation; the constraints provided by astronomical observations 
are very indirect. It is prudent to put these theories to a laboratory test. 

Secondly, computer speed limits the range of spatial scales that can be modeled in the sim- 
ulations. Barring unforseen algorithmic breakthroughs, the smallest resolvable scale in a three- 
dimensional simulation improves only as the fourth root of the rate of arithmetic operations. Here 
it must be acknowledged that the large magnetic diffusivity of liquid metals severely limits the num- 
ber of degrees of freedom in the magnetic field that can be excited. The simulations are well ahead 
of any forseeable experiment in this respect. In fact, simulations indicate that magnetic Reynolds 
numbers and Lundquist numbers at least 100 times larger than the minimum necessary for linear 
instability are required for dynamo action in the absence of an externally imposed field parallel 
to the rotation axis (Fleming, Stone, & Hawley 2000). On the other hand, the viscous Reynolds 
number of such an experiment would be R c > 10 6 , a value still out of reach of direct numerical 
simulations. Also, small R m need not restrict the experiment to linear behaviors. In the local disk 
simulations of Fleming, Stone, & Hawley (2000) at about twice the minimum R m for linear MRI, a 
violently fluctuating nonlinear state was reached in which the time-averaged magnetic energy was 
about 25 times larger than that of the externally imposed field. 

Although large R m is the rule in astrophysics, the dimensionless parameters of some systems 
may be similar to those of our proposed experiment, viz. R m and S of order unity, R e very large, 
and an externally imposed field. Such systems include the inner parts of relatively cool disks 
(protostellar disks and quiescent cataclysmic variables, for example) around stars with their own 
magnetic moments (Gammie 1996; Gammie &: Menou 1998). 
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Lastly, relatively little laboratory MHD work has been done in which the inertia of the fluid 
is important (large plasma (3). The experimental field appears somewhat underdeveloped when 
measured against its potential importance to geophysics and astrophysics. Because it promises to 
be achievable at fairly modest cost in a classic experimental framework (Couette flow), MRI is a 
good place to start. 

This work was supported by the U.S. Department of Energy [H.J.] and by NASA grant NAG5- 
8385 [J.G.] 



A. Derivation of linearized equations 

For completeness, eqs. (2)-(6) are derived here, although much the same derivation can be 
found in Chandrasekhar (1961). The equations of incompressible MHD are 

B + v VB-B Vv = r/V 2 B, VB = 0, 

! B VB o 

v + v-Vv + p L VP = zA7 2 v, V-v = 0, 

Airp 

in which P = p + B 2 /87r, is the hydrodynamic plus magnetic pressure. In cylindrical coordinates, 
near an equilibrium B = Be z = constant and v = rQ(r)eQ, linearized axisymmetric perturbations 
<5v and <5B satisfy 

6B r -Bd z 6v r = ri(d r d} + d 2 z )5B r ., (Al) 
5B e - Bd z 5v e - 5B r rd r n = r)(d r dl + d 2 )5B e , (A2) 

ID 73 

5v r -2M6v e + dr - A — d z 5B r = v{d r d\ + d 2 z )5v r , (A3) 

p 4lTp 

5v e + 5v r dt{rVL) - -^-d z 5B e = v(d r dl + d 2 z )5v , (A4) 

47T/9 

5v z +d z —-^-8 z 5B z = u(dld r + d 2 )8v z , (A5) 
p <±7rp 

d}5B r + d z 5B z = 0, (A6) 
dl&v r + d z 8v z = 0, (A7) 

in which the dot denotes d/dt, and other recurring operators are 

d „ d 



d z = —, d r = — , d} = — +-. 
oz or or r 



Eqs. (A3) & (A5) presume that p, like rj and v, is spatially constant. Applying d\ to eq. (A3) and 
d z to eq. (A5) and summing the results, one finds that 

xp 

(dld r + d 2 )^ = dl(2n5v e ), 
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in view of eqs. (A6) & (A7). With another application of d r , this becomes 

SP 

(d r dl + dl)U = d^(2nSv e ), where II = 2Sl5v e - d r — , (A8) 

so that the radial euler equation (A3) can be stated as 

5v r -U- — d z 5B r = v(d r dl + d 2 z )8v r . (A9) 

With the z dependences given by eq. (1) for the linearized quantities, eqs. (A2), (A4), (Al), (A9), 
and (A8) reduce to eqs. (2), (3), (4), (5), and (6), respectively. 



REFERENCES 

Balbus, S. A. & Hawley, J. F. 1991 Astrophys. J. 376, 214. 
Balbus, S. A. & Hawley, J. F. 1992 Astrophys. J. 392, 662. 
Brahme, A. 1970 Physica Scripta 2, 108. 
Chandrasekhar, S. 1960 Proc. Nat. Acad. Sci. 46, 253. 

Chandrasekhar, S. 1961 Hydrodynamic and Hydro-magnetic Stability. (London: Oxford University 
Press). 

Chang, T. S. & Sartory, W. K. 1967 Proc. R. Soc. Lond. A 301, 451. 
Chen, C.-K. & Chang, M. H. 1998 J. Fluid Mech. 366, 135. 
Donnelly, R. J. & Caldwell, D. R. 1964 J. Fluid Mech. 19, 257. 
Donnelly, R. J. & Ozima, M. 1960 Phys. Rev. Lett. 4 (10), 497. 
Donnelly, R. J. & Ozima, M. 1960 Proc. R. Soc. Lond. A 266, 272. 
Fleming, T.P., Stone, J.M., Hawley, J.F. ApJ, 530, 464. 
Gammie, C. F. 1996 ApJ, 457, 356. 
Gammie, C. F. & Menou, K. 1998 ApJ, 492, L75. 

Hassard, B. D., Chang, T. S. & Ludford, G. S. S. 1972 Proc. R. Soc. London A 327, 269. 
Ji, H., Goodman, J., & Kageyama, A. 2001 astro-ph/0103226 [Paper I]. 
Terquem, C. k Papaloizou, J. C. B. 1996 Mon. Not. R. Astr. Soc. 279, 767. 
Soundalgekar, V. M., Ali, M. A., & Takhar, H. S. 1994 Int. J. Energy Res. 18, 689. 



-14- 



Takhar, H. S., Ali, M. A., & Soundalgekar, V. M. 1994 Appl. Sci. Res. 48, 1. 

Vislovich, A. N., Novikov, V. A., & Sinitsyn, V. A. 1986 J. Appl. Mech. Tech. Phys. 27 (1), 72. 

Velikhov, E. P. 1959 Sov. Phys. JETP 36 (2), no. 5, 995. 



This preprint was prepared with the AAS IATjrjX macros v5.0. 



-15- 




Fig. 1. — Marginal stability for liquid gallium Couette flow between conducting cylinders of radii 
r\ = 5 cm, r-2 = 15cm, and height h = 10cm. Solid lines computed from eqs. (2)-(6); dashed from 
inviscid approximation (26)- (27). Lower curves for dimensionless vorticity ( = 2/11, upper ones 
for ( = 4/7. Instability occurs above the curves. In dimensionless parameters, S « 0.92(5/10 kG); 
and R m « 0.66(ft 2 /100 rad s" 1 ) (upper curves), R m « 0.73(^ 2 /100rad s" 1 ) (lower). 
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Fig. 2. — Like Fig. (1), but for insulating cylinders (8). Viscous and inviscid results differ by less 
than the line thickness. 
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Fig. 3. — Global and closeup views of marginal eigenmode with conducting boundaries. B = 3kG, 
Sli = 314., ft 2 = 37.9 rad s" 1 & ( = 0.0632. Solid curve: f3 r . Short-dashed: (3 Z Dot-long-dashed: 
(3q x 5. Dotted: ip r x 1/3. Long-dashed: Dot-short-dashed: ip z x 0.07. 
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Fig. 4. — Like Fig. 3 but for insulating boundaries and Qi = 284., 0, 2 = 34.4 rad s _1 , C = 0.0632. 
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Fig. 5. — Narrow-gap modes for P m = and e = 1. Left panel: Curves of marginal stability for 
C = 2/11 (lower curves) and £ = 4/7 (upper). Solid line for conducting walls, dashed for insulating, 
and dotted for local approximation (28). Right panel: Narrow-gap eigenfunctions (3 r (solid curves) 
and /3g (dashed), for £ = 2/11 and S = 0.4. Since eigenfunctions are symmetric about center of 
gap (x = 0), only half of each is shown: conducting on left, insulating on right. 
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Fig. 6. — Visualizations of the MRI eigenmodes for the Rayleigh-stable cases from Table 1 at 
B Zj o = 3 kG. Left: conducting boundaries. Right insulating. Solid and dotted lines indicate 
positive and negative values, respectively. See eq. (31) for definitions of flux and stream functions 
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Fig. 7. — Rayleigh-unstable cases from Table 1 at of B Zi q 
conducting boundaries. 



= kG (left) and 3 kG (right), both with 
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Table 1. Growth rates in gallium 



fii = 413.6, ft 2 = 50. rad s" 1 
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